Accelerated dynamics with the dynamical activation-relaxation technique 
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The dynamics of many atomic systems is controlled by activated events taking place on a time 
scale which is long compared to that associated with thermal vibrations. This often places problems 
of interest outside the range of standard simulation methods such as molecular dynamics. We 
present here an algorithm, the dynamical activation-relaxation technique (DART), which slows down 
thermal vibrations, while leaving untouched the activated processes which constitute the long-time 
dynamics. As an example, we show that it is possible to accelerate considerably the dynamics of 
self-defects in a 1000-atom cell of c-Si over a wide range of temperatures. 
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Developing algorithms that stretch the time scale ac- 
cessible to computer simulations has been a major chal- 
lenge in computational physics, chemistry and biology. 
Standard techniques such as molecular dynamics are 
strongly constrained by the presence of high-frequency 
modes in dense materials and must therefore use an in- 
tegration time step on the order of a femtosecond, limit- 
ing the simulation length to microseconds at best. Over 
the years, a number of accelerated algorithms have been 
proposed for discrete systems where it is possible to enu- 
merate, in advance, all the possible barriers. The original 
ideas for these algorithms are due to Bortz, Kalos and 
Lebowitz |l|. They were first applied to MBE growth in 
1986 by Voter 0. This method and many variants are 
now common simulation techniques in surface science. 

Algorithmic progress has been much slower for sys- 
tems in which either the number of pathways, or their 
activation energies, change all the time. Many meth- 
ods have been proposed for sampling events efficiently 
H II II 11 Si, IE E3 but onl y a few of these can provide 
a time scale: hyper- molecular dynamics and related ap- 
proches 0, 0] as well as temperature-assisted dynam- 
ics, introduced by Sorensen and Voter are limited 
to systems with a relatively narrow distribution of bar- 
riers and cannot be used easily at high temperatures or 
on generic problems. Other methods, such as a kinetic 
Monte Carlo scheme with on-the-fly calculation of the 
barriers [T^j generate impressive acceleration. However, 
their application is limited to relatively simple systems, 
and entropic effects are not fully included, hence detailed 
balance cannot be ensured. 

In this paper, we present an algorithm, the dynami- 
cal activation-relaxation technique (DART), with a self- 
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correcting accelerating factor that can reach a few orders 
of magnitude at physically relevant temperatures while 
overcoming some of the limitations of these accelerated 
methods. It is based on the thermodynamically-weighted 
activation-relaxation technique (THWART) which we in- 
troduced recently , and combines molecular dynamics 
with Monte Carlo methods to reach a dynamically correct 
acceleration of the slow processes in complex materials. 
In order to assess the efficiency of DART, we apply the 
method to the diffusion of vacancies and interstitials in 
c-Si. 

The energy landscape of systems with dynamics con- 
trolled by rare events can be divided into two types of 
subregions, as shown in Fig.^ the basins, closed regions 
of the energy landscape, in which the system is confined 
for extended periods and which can contain many local 
minima, and the "activated part of phase space", sam- 
pled only when a rare fluctuation pushes the system from 
one basin to another. Following THWART, we delineate 
the basins based on the value of the lowest local curva- 
ture of the energy landscape (i.e., the lowest eigenvalue 
of the hessian matrix); any configuration with a lowest 
curvature below the threshold is considered to be in the 
saddle region. The exact threshold value depends on the 
system studied as well as on the simulation temperature. 

The simulation starts with standard molecular dynam- 
ics in a microcanonical ensemble. After equilibration, 
the MD run is pursued but with regular evaluation of 
the lowest local curvature Ao- As long as Ao remains 
above a threshold value At, the system is considered to 
reside in the original basin. As soon as this threshold is 
reached, however, the activation phase of THWART is 
launched: the MD simulation is stopped, the current ve- 
locities stored and the system is moved from one basin to 
another. The step is designed to bring the configuration 
from the edge of a basin to that of a neighboring one, in 
a fully reversible manner and at constant potential en- 
ergy in order to respect detailed balance; it is identical to 



2 




FIG. 1: The energy landscape of systems with a long-time 
dynamics determined by activated events. 



that used in THWART. Each activation move consists of 
a sequence of small steps with size Ax, defined as follows: 

Xi+i =Si + ^- (hi+i + hi) + c (P^ + F±) , (1) 

where h is the direction of lowest curvature and F is 
the force perpendicular to this vector; c is a scalar tuned 
at each step to keep the total configurational energy con- 
stant. The orientation of ho is chosen away from the 
basin, while successive orientations are chosen such that 
hi + i - hi > 0. This step is repeated until the lowest curva- 
ture of the energy landscape reaches again the threshold 
value A t . Note that each individual step, and hence the 
whole trajectory, is fully reversible. Since each activation 
move connects two points in phase space with the same 
potential energy and since the Jacobian of transformation 
is equal to one detailed balance is fully respected. 
Within THWART, activated moves are always accepted; 
the MD is therefore continued from the end of the acti- 
vation path using the stored velocities, ready for a new 
activation. 

THWART ensures a proper thermodynamical sam- 
pling of the phase space; the THWART trajectory, how- 
ever, does not follow the real dynamics. In particular, 
it crosses highly activated pathways as easily as those 
with a low activation barrier. In order to recover the 
actual dynamics, it is necessary to first determine the 
height of the barriers crossed. By construction, activa- 
tion pathways have constant energy and the deformation 
energy stored in the degrees of freedom not directly in- 
volved in the activation process is used as a bath to en- 
force this constraint. It is nevertheless straightforward 
to to determine the energy associated with the diffusion 
path: the total force is split into a component parallel to 
the direction of lowest curvature (F^) and a component 
perpendicular to this direction (F^), at each step along 
the activation trajectory; using the first component, the 
change in total energy due to the parallel displacement 



is written 

AE^ =-fm -dR. (2) 

The activation energy barrier (from the edge of the basin) 
is then defined as the change in energy from the edge of 
the basin to the maximum energy change projected along 
this path. 

The probability that the kinetic energy at the begin- 
ning of the activation trajectory suffices to bring the con- 
figuration over the barrier through this highest point is 
given by exp(— /3max(Ai?")), with inverse temperature 
(3 = l/(kbT). Assuming around each first-order saddle 
point a quadratic behavior of the energy in the transition 
plane, the energy at the nearest saddle point will on av- 
erage be \k b T lower. We therefore take for the activation 
barrier, faced from the basin boundary: 

A£ act = max(A£" ) - h b T. (3) 

To retrieve the dynamics in a statistically correct man- 
ner, we should accept the activation move with probabil- 
ity exp(— /3A_E act ), and otherwise continue in the original 
basin. Typically in a system with activated dynamics, 
these acceptance probabilities are rather small, and the 
system will bounce back and forth in the basin many 
times before eventually escaping from it; hence the slow 
dynamics. However, the speed of the simulation can be 
enhanced by a (constant) nominal boost factor X b reduc- 
ing the number of such bounces. The acceptance proba- 
bility for an activation pathway then becomes 

P cmss = min [l,X b exp (-/3A£ act )] . (4) 

If, on top of boosting the acceptance probabilities by 
a factor of X b , the time scale is stretched by the same 
factor, the long-time dynamics is untouched, provided 
it is indeed determined by activated processes with bar- 
riers exceeding \m(Xf,)kbT, while the suppression of the 
in-basin dynamics (by a factor X b ) as well as the suppres- 
sion of less activated processes (by a smaller factor) does 
not affect the long-time dynamics. Once barriers below 
ln(X b )k b T are encountered, inevitably some distortion of 
the dynamics occurs. To alleviate the distortion some- 
what, if we encounter such a low barrier, we stretch the 
time scale since the previous event by an on-the-fly cor- 
rected boost factor 

X c g = Pcross • exp(j3AE act ), (5) 

rather than the nominal boost factor X b ; this recovers 
correct dynamics for systems in which the activation en- 
ergy is constant, even if the chosen nominal boost factor 
is too large. 

A DART simulation proceeds in the following se- 
quence: (1) at time t = 0, a microcanonical molecular 
dynamics simulation is launched and the value of the 
lowest local curvature Ao is monitored at regular inter- 
vals (typically, every 50 steps); (2) when Aq reaches the 
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threshold value At, the molecular dynamics is stopped 
and the velocities arc saved; (3) following Eq. itera- 
tively, the configuration is brought into a new basin and 
the activation energy A.E ac t is computed along the ac- 
tivated pathway; (4) the event is accepted with prob- 
ability Pcross) (5) if the event is accepted, the time is 
incremented by At = £md ■ A g, in which £md is the time 
spent doing molecular dynamics since the previous ac- 
cepted event, and the molecular dynamics simulation is 
continued starting at the new edge with the same veloc- 
ities. If it is rejected, the molecular dynamics simulation 
is simply continued from the initial edge. 

To demonstrate the efficiency of DART, we consider 
the diffusion of vacancies and interstitials in c-Si, de- 
scribed by the Stillinger- Weber potential. Both tvpes of 
defects have been well characterized previously [3 EH • 
Vacancy diffusion is associated with a single activation 
barrier of 0.43 eV 0|. An interstitial can take four dif- 
ferent stable topologies: tetrahedral, hexagonal, bond- 
centered and split or dumbbell. It diffuses through many 
mechanisms with activation barriers between 0.65 and 
1.62 eV [16|. Therefore, these two defects provide us 
with various levels of complexity to test DART. 

Figure [^shows an Arrhenius plot of the diffusion rates 
obtained by MD and by DART with various boost fac- 
tors. We characterize the diffusion rates by the hopping 
rates of the defects, which gives better statistics than the 
mean squared displacements per unit of time. As can be 
seen, an excellent agreement between the two methods 
is achieved for both the vacancy, characterized by a sin- 
gle energy barrier, and for the interstitial, which shows a 
more complex behavior. 

As mentioned above, DART adjusts automatically for 
barriers lower than ln(Af,)fcf,T. Table Ogives the nom- 
inal and effective boost factors for the various DART 
simulations plotted in Fig. [21 As expected, the gain 
in efficiency increases rapidly with decreasing tempera- 
ture. From a factor of 1 at about HOOK (not shown), 
the effective boost for interstitial diffusion reaches 10 at 
900 K, 74 at 750 K and almost 150 at 600 K, an experi- 
mentally relevant temperature. The extra computational 
effort in DART is largely due to the computation of the 
lowest curvature, both at regular intervals and during the 
THWART events. Averaged over the whole simulation, 
with an evaluation at every 50 MD steps, a DART time 
step costs slightly less than two MD time steps. 

In conclusion, we have presented here an accelerated 
molecular dynamical method the dynamical activation- 
relaxation method (DART). This algorithm provides a 
tunable acceleration parameter that can be adjusted to 
suit the specific problems studied. In addition to provid- 
ing a significant acceleration over MD, DART has numer- 
ous advantages: (1) the algorithm is not very sensitive to 
the various parameters - it can automatically correct for 
boost factors an order of magnitude or more too large; (2) 
it computes relaxation trajectories and activation barri- 
ers on the fly, leading to a very low overhead (on average, 
a DART time step is less than twice the cost of an MD 
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FIG. 2: Diffusion rate as a function of inverse temperature 
for vacancy (top) and the interstitial (bottom) diffusion in Si. 
The open circles represent MD results and the squares DART 
results. Error bars are obtained from the square root of the 
total number of events. The nominal and averaged effective 
boost factors for these simulations are given in Table [I] 



TABLE I: Details of the various DART simulations 



Defect 


Temperature 


A 


Boost 




(K) 


(eV/ A 2 ) Nominal Effective 


Vacancy 


900 


-5 


6 


3.6 




900 


-5 


30 


6.7 




600 


-5 


6 


5.4 




600 


-5 


30 


19 




450 


-5 


600 


271 


Interstitial 


900 


-7 


60 


10 




750 


-7 


600 


74 




600 


-7 


1200 


148 
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time step); (3) DART is not slowed down by the local 
rearrangements which take place in the basin (i.e, below 
threshold) — contrary to other accelerated methods, it is 
therefore possible to use DART to accelerate the dynam- 
ics of more complex systems such as glasses and proteins; 
(4) since the events are easily labeled, it is possible to use 
various tricks to avoid repeating the same event over and 
over again - this can lead to a large increase in efficiency; 
finally, (5) the limits of DART are well behaved: a zero- 
boost reduces DART to standard MD, while an infinite 
boost factor recovers THWART and still ensures proper 
thermodynamical sampling. 



Tests on a vacancy and an interstitial in c-Si have 
shown that this method remains accurate dynamically 
with an increased computational efficiency of two or more 
orders of magnitude, at temperatures as high as 450 K for 
the vacancies or 600 K for interstitial diffusion. As men- 
tioned above, however, this method can also be applied to 
much more complex situations where current accelerated 
algorithms fail. 
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